{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic回归作业 利用Logistic回归技术实现糖尿病发病预测\n",
    "\n",
    "## 一、\t数据说明： Pima Indians Diabetes Data Set（皮马印第安人糖尿病数据集） 根据现有的医疗信息预测5年内皮马印第安人糖尿病发作的概率。\n",
    "数据链接：https://archive.ics.uci.edu/ml/datasets/Pima+Indians+Diabetes\n",
    "p.s.: Kaggle也有一个Practice Fusion Diabetes Classification任务，可以试试:)\n",
    "https://www.kaggle.com/c/pf2012-diabetes\n",
    "\n",
    "1)\t文件说明\n",
    "pima-indians-diabetes.csv：数据文件\n",
    "\n",
    "2)\t字段说明\n",
    "数据集共9个字段: \n",
    "pregnants：怀孕次数\n",
    "Plasma_glucose_concentration：口服葡萄糖耐量试验中2小时后的血浆葡萄糖浓度\n",
    "blood_pressure：舒张压，单位:mm Hg\n",
    "Triceps_skin_fold_thickness：三头肌皮褶厚度，单位：mm\n",
    "serum_insulin：餐后血清胰岛素，单位:mm\n",
    "BMI：体重指数（体重（公斤）/ 身高（米）^2）\n",
    "Diabetes_pedigree_function：糖尿病家系作用\n",
    "Age：年龄\n",
    "Target：标签， 0表示不发病，1表示发病\n",
    "\n",
    "## 二、作业要求：\n",
    "1.\t对数据做数据探索分析（可参考0_EDA_ diabetes.ipynb，不计分）\n",
    "2.\t适当的特征工程（可参考1_FE_ diabetes.ipynb，不计分）\n",
    "3.\t采用5折交叉验证，分别用log似然损失和正确率，对Logistic回归模型的正则超参数调优。（各50分）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pregnants</th>\n",
       "      <th>Plasma_glucose_concentration</th>\n",
       "      <th>blood_pressure</th>\n",
       "      <th>Triceps_skin_fold_thickness</th>\n",
       "      <th>serum_insulin</th>\n",
       "      <th>BMI</th>\n",
       "      <th>Diabetes_pedigree_function</th>\n",
       "      <th>Age</th>\n",
       "      <th>Target</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.639947</td>\n",
       "      <td>0.866045</td>\n",
       "      <td>-0.031990</td>\n",
       "      <td>0.670643</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>0.166619</td>\n",
       "      <td>0.468492</td>\n",
       "      <td>1.425995</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.844885</td>\n",
       "      <td>-1.205066</td>\n",
       "      <td>-0.528319</td>\n",
       "      <td>-0.012301</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>-0.852200</td>\n",
       "      <td>-0.365061</td>\n",
       "      <td>-0.190672</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.233880</td>\n",
       "      <td>2.016662</td>\n",
       "      <td>-0.693761</td>\n",
       "      <td>-0.012301</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>-1.332500</td>\n",
       "      <td>0.604397</td>\n",
       "      <td>-0.105584</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.844885</td>\n",
       "      <td>-1.073567</td>\n",
       "      <td>-0.528319</td>\n",
       "      <td>-0.695245</td>\n",
       "      <td>-0.540642</td>\n",
       "      <td>-0.633881</td>\n",
       "      <td>-0.920763</td>\n",
       "      <td>-1.041549</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.141852</td>\n",
       "      <td>0.504422</td>\n",
       "      <td>-2.679076</td>\n",
       "      <td>0.670643</td>\n",
       "      <td>0.316566</td>\n",
       "      <td>1.549303</td>\n",
       "      <td>5.484909</td>\n",
       "      <td>-0.020496</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   pregnants  Plasma_glucose_concentration  blood_pressure  \\\n",
       "0   0.639947                      0.866045       -0.031990   \n",
       "1  -0.844885                     -1.205066       -0.528319   \n",
       "2   1.233880                      2.016662       -0.693761   \n",
       "3  -0.844885                     -1.073567       -0.528319   \n",
       "4  -1.141852                      0.504422       -2.679076   \n",
       "\n",
       "   Triceps_skin_fold_thickness  serum_insulin       BMI  \\\n",
       "0                     0.670643      -0.181541  0.166619   \n",
       "1                    -0.012301      -0.181541 -0.852200   \n",
       "2                    -0.012301      -0.181541 -1.332500   \n",
       "3                    -0.695245      -0.540642 -0.633881   \n",
       "4                     0.670643       0.316566  1.549303   \n",
       "\n",
       "   Diabetes_pedigree_function       Age  Target  \n",
       "0                    0.468492  1.425995       1  \n",
       "1                   -0.365061 -0.190672       0  \n",
       "2                    0.604397 -0.105584       1  \n",
       "3                   -0.920763 -1.041549       0  \n",
       "4                    5.484909 -0.020496       1  "
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "%matplotlib inline\n",
    "\n",
    "train = pd.read_csv('FE_pima-indians-diabetes.csv')\n",
    "train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 准备数据\n",
    "y_train = train['Target']\n",
    "X_train = train.drop([\"Target\"], axis = 1)\n",
    "\n",
    "# 保存特征名字以备后用（可视化）\n",
    "feat_name = X_train.columns\n",
    "\n",
    "# sklearn的学习器大多之一稀疏数据输入，模型训练会快很多\n",
    "# 查看一个学习期是否支持稀疏数据，可以看fit函数是否支持：X:{array-like, sparse matrix}\n",
    "# 可自行用timeit比较稠密数据和系数数据的训练时间\n",
    "from scipy.sparse import csr_matrix\n",
    "x_train = csr_matrix(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logloss of each fold is : [0.48797856 0.53011593 0.4562292  0.422546   0.48392885]\n",
      "cv logloss is: 0.47615970944434044\n"
     ]
    }
   ],
   "source": [
    "# 使用默认参数的Logistic Regression\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "lr = LogisticRegression(solver='liblinear')\n",
    "\n",
    "# 使用交叉验证用于评估模型性能和进行参数调优（模型选择）\n",
    "# 分类人物中交叉验证缺省是采用 StratifiedKFold\n",
    "# 数据集比较大，采用5折交叉验证\n",
    "from sklearn.model_selection import cross_val_score\n",
    "loss = cross_val_score(lr, X_train, y_train, cv=5, scoring = 'neg_log_loss')\n",
    "print('logloss of each fold is :', -loss)\n",
    "print('cv logloss is:', -loss.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 正则化的logistic Regression 及参数调优\n",
    "logistic回归的需要调整超参数有：C（正则系数，一般在log域（取log后的值）均匀设置候选参数）和正则函数penalty（L2/L1） \n",
    "目标函数为：J =  C* sum(logloss(f(xi), yi)) +* penalty \n",
    "\n",
    "在sklearn框架下，不同学习器的参数调整步骤相同：\n",
    "设置参数搜索范围\n",
    "生成GridSearchCV的实例（参数）\n",
    "调用GridSearchCV的fit方法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=5, error_score='raise-deprecating',\n",
       "             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,\n",
       "                                          fit_intercept=True,\n",
       "                                          intercept_scaling=1, l1_ratio=None,\n",
       "                                          max_iter=100, multi_class='warn',\n",
       "                                          n_jobs=None, penalty='l2',\n",
       "                                          random_state=None, solver='liblinear',\n",
       "                                          tol=0.0001, verbose=0,\n",
       "                                          warm_start=False),\n",
       "             iid='warn', n_jobs=None,\n",
       "             param_grid={'C': [0.001, 0.01, 0.1, 1, 10, 100, 1000],\n",
       "                         'penalty': ['l1', 'l2']},\n",
       "             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,\n",
       "             scoring='accuracy', verbose=0)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "#需要调优的参数\n",
    "# 请尝试将L1正则和L2正则分开，并配合合适的优化求解算法（slover）\n",
    "#tuned_parameters = {'penalty':['l1','l2'],\n",
    "#                   'C': [0.001, 0.01, 0.1, 1, 10, 100, 1000]\n",
    "#                   }\n",
    "penaltys = ['l1','l2']\n",
    "Cs = [0.001, 0.01, 0.1, 1, 10, 100, 1000]\n",
    "tuned_parameters = dict(penalty = penaltys, C = Cs)\n",
    "\n",
    "lr_penalty= LogisticRegression(solver='liblinear')\n",
    "grid_neg_log_loss= GridSearchCV(lr_penalty, tuned_parameters,cv=5, scoring='neg_log_loss')\n",
    "grid_accuracy = GridSearchCV(lr_penalty, tuned_parameters, cv=5, scoring='accuracy')\n",
    "grid_neg_log_loss.fit(X_train,y_train)\n",
    "grid_accuracy.fit(X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4760244986295943\n",
      "{'C': 1, 'penalty': 'l1'}\n",
      "0.7747395833333334\n",
      "{'C': 0.1, 'penalty': 'l2'}\n"
     ]
    }
   ],
   "source": [
    "# examine the best model\n",
    "print(-grid_neg_log_loss.best_score_)\n",
    "print(grid_neg_log_loss.best_params_)\n",
    "print(grid_accuracy.best_score_)\n",
    "print(grid_accuracy.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEICAYAAABF82P+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3hUVfrA8e87SUhCEmoiLSBFegu9iFhQYEVBBZbiCoiKrLK67upadtcuq6srqz9dFSkiKkWxoKKgWBCVTlCagEiJoECoAZMQ5v39cSchhAkZMjOZTPJ+nmeezL333DvvJeG+c+459xxRVYwxxpiCXKEOwBhjTOlkCcIYY4xXliCMMcZ4ZQnCGGOMV5YgjDHGeGUJwhhjjFchSRAi0ldEfhCRLSJyj5ftE0Qk1fPaJCIHQxGnMcaUZ1LSz0GISASwCbgMSAOWA8NUdX0h5f8EtFPV0Wc6bmJiotavXz/A0RpjTNm2cuXKfaqa5G1bZEkHA3QGtqjqVgARmQkMALwmCGAY8EBRB61fvz4rVqwIWJDGGFMeiMj2wraF4hZTHWBnvuU0z7rTiMi5QAPgs0K2jxGRFSKyYu/evQEP1BhjyrNQJAjxsq6w+1xDgbdU9YS3jao6UVU7qmrHpCSvNSSfZOe4i72vMcaUVaFIEGlA3XzLycCuQsoOBWYEM5jt6Ufp9NinPPDeWn745UgwP8oYY8JKKNoglgONRaQB8DNOEhhesJCINAWqAt8GMxhVuLhpEjOW72Tat9tpX68Kw7ucS7/WtYitEBHMjzYmJI4fP05aWhqZmZmhDsWUoJiYGJKTk4mKivJ5nxLvxQQgIpcD/wUigCmq+piIPAysUNW5njIPAjGqelo3WG86duyo/jRSHziazZxVabyxbAdb9x6lUkwk17RPZljnejStmVDs4xpT2vz0008kJCRQvXp1RLzd8TVljaqSnp7OkSNHaNCgwSnbRGSlqnb0tl9IEkQw+JsgcqkqS3/az4xlO/jo+1/IPuGmw7lVGda5Hle0qUVMlNUqTHjbsGEDzZo1s+RQzqgqGzdupHnz5qesP1OCCMUtplJNROjasDpdG1bngSuzedtTq7jzzTU8/P46q1WYMsGSQ/lTnN+5JYgzqBZXgRsvaMgNPRrk1SreWLqDV77ZZrUKU64MeclpCpx1c7cQR2JKko3F5IPcWsUzQ9ux5L5e/KNfcw4cy+bON9fQ+bFPeXDuOjb9aj2gjPFVfHx83vu+fftSpUoVrrjiikLLjxo1ijp16pCVlQXAvn378GfkhFGjRvHWW28Va99XXnmFcePGFfuzC9OlSxdSUlKoV68eSUlJpKSkkJKSwrZt287qOG+//TYbN24MSExWgzhLBWsVbyw9tVYxvHM9+lmtwhif3XXXXRw7doyXXnrpjOUiIiKYMmUKf/zjH0sospK1dOlSwElAK1as4LnnnivWcd5++21cLhfNmjXzOyarQRRTbq3i2WFOreLvlzfnwNFs/mq1CmPOSq9evUhIKLpN789//jMTJkwgJyfntG1PPvkknTp1ok2bNjzwwMmReR555BGaNWvGZZddxrBhw3jqqadO23fhwoW0a9eO1q1bM3r06Lxayrx582jWrBk9evTgtttu81rD2b59O7169aJNmzb06tWLHTt2APDmm2/SqlUr2rZtS8+ePQFYt24dnTt3JiUlhTZt2rB582bf/oGAjz76iG7dutG+fXuGDBnC0aNHASe5tmjRgjZt2nD33Xfz1VdfMW/ePO64445i1T4KshpEAFSLq8BNPRty4wVWqzDh5aH317F+1+Eiy63f7ZTJbYs4kxa1K/HAlS39jq2gevXq0aNHD6ZPn86VV16Zt37BggVs3ryZZcuWoar079+fRYsWUbFiRebMmcPq1avJycmhffv2dOjQ4ZRjZmZmMmrUKBYuXEiTJk0YMWIEL7zwAmPHjuXmm29m0aJFNGjQgGHDhnmNady4cYwYMYKRI0cyZcoUbrvtNt59910efvhh5s+fT506dTh40BmM+sUXX+T222/n2muvJTs7mxMnvA4QcZo9e/bw+OOPs3DhQipWrMhjjz3GM888ww033MC8efNYt24dIsLBgwepUqUKl19+OYMGDeKqq64q5r/0SVaDCCCrVRgTXPfddx9PPvkkbvfJ4XEWLFjAggULaNeuHe3bt2fjxo1s3ryZxYsXM2DAAGJjY0lISDglqeT64YcfaNCgAU2aNAFg5MiRLFq0iI0bN9KwYcO8ZwYKSxDffvstw4c7z/led911LF68GIDzzz+fUaNG8fLLL+clgm7dujF+/HieeOIJtm/fTmxsrE/n/M0337B+/Xq6d+9OSkoKr7/+Otu2baNatWq4XC5uuukm3nnnHeLi4nz8V/Sd1SCCJH+tYsnWU3tAdfT0gLJahQk1X7/pl5ZeTOeddx4pKSnMnj07b52qcu+993LzzTefUnbChAlFHq+w58CK+3xYblfSF198kaVLl/Lhhx+SkpJCamoqw4cPp0uXLnz44Yf06dOHSZMmcckll/gUY9++fZk+ffpp21asWMEnn3zCzJkzeeGFF1iwYEGx4i6M1SCCTETo1ujUWsV+q1UYU2x///vfT2lL6NOnD1OmTCEjIwOAn3/+mT179tCjRw/ef/99MjMzycjI4MMPPzztWM2aNWPbtm1s2bIFgOnTp3PhhRfSrFkztm7dmncPf9asWV5j6d69OzNnzgTg9ddfp0ePHgD8+OOPdOnShYcffpjExER27tzJ1q1badiwIbfddhv9+/fnu+++8+l8u3fvzpdffsnWrVsBOHr0KJs3b+bIkSMcPnyYK664ggkTJrB69WoAEhISOHIkMNcUq0GUoIK1ijcK1CqGd6nH5a2tVmHKlwsuuICNGzeSkZFBcnIykydPpk+fPoWWb9myJe3bt2fVqlUA9O7dmw0bNtCtm1O7iY+P57XXXqNTp07079+ftm3bcu6559KxY0cqV658yrFiYmKYOnUqgwcPJicnh06dOjF27Fiio6P53//+R9++fUlMTKRz585eY3n22WcZPXo0Tz75JElJSUydOhVwGo83b96MqtKrVy/atm3L448/zmuvvUZUVBQ1a9bk/vvv9+nfp0aNGkyePJkhQ4aQnZ0NwPjx44mNjeWaa64hKysLt9vN008/DTi3w26++Wb+85//8O677/rVHdiG2gix9Iws5qxKY8aynfy07+QYUNd2qUfjGva0tgm8DRs2nDbcQlFKyy2ms5WRkUF8fDzHjh2jZ8+eTJw4kfbt25/VvqrKrbfeSuPGjbnjjjuCHHFwefvd21AbpVj1+GjG9GzETRc05Nut6cxYtpPXl27nlW+20am+01ZhtQoTauGWGHKNGTOG9evXk5mZyciRI31ODgAvv/wy06ZNIzs7m3bt2p3WxlEeWA2iFCpYq6gcG8U17eswvLPVKoz/ilODMGWD1SDKAG+1iteWbGfq11arMMaUHEsQpZiI0L1RIt0bJZKe0SKvVvGX2Wt46P31VqswxgSVJYgwUbBW8cbSHafVKmYs3YHLJWF7v9gYU7pYgggzp9YqsnhrZRozlu3gL7PXEOESaiREc8KtRLhsvH8TQFP7OT+vP/1ZAlN22YNyYax6fDQ3X9iIz++8iDdu6kLlmEh2Hcrkn++tLfaToMaUhNzhvlNTU+nWrRstW7akTZs2hT6QZsN9+86G+w60MP92lFuraFwjgR37j/HG0h0kxUdzx2VNQh2aMWdUsWJFXn31VRo3bsyuXbvo0KEDffr0oUqVKqeVteG+fWPDfQfSsf2wZx0cPxbqSPw26+ZufHPPJfy+YzLPLNzM9CXbQx2SMWfUpEkTGjduDEDt2rU555xz2Lt3r9eyNty3Dfdd8g5sg6wjsHsNbPgAmhc+q1U4EBHGX92a9Ixs7n9vLdXjKnB561qhDsuUVh/dA798X3S5XzzjBuXWts+kZmv43eNnHcqyZcvIzs6mUaNGXrfbcN823HfJq9MeaqVAVCzMuhYWPgJu335xpVVkhIvnhrenQ72q/HlmKt/8uC/UIRlzRrt37+a6665j6tSpuFyFX5ZsuG8b7rvk3bQQjmfCvL/CV0/BrtUwcBJUrBbqyIottkIEk0Z25PcvfcuYV1cyc0xXWtWpXPSOpnzx9Zt+ENvpDh8+TL9+/Xj00Ufp2rXrGcvacN+nsuG+S0pUDPR/Dq74L/y0CCZeBLt9G463tKpSsQLTRnemUkwko6YuZ0d6+LezmLIlOzubq6++mhEjRjB48GCf9rHhvktuuG9LEPmJQMfr4fqP4MRxmNwb1nj/wwgXtSrH8uoNnclxu7luylL2HskKdUjG5Jk9ezaLFi3ilVdeyevWmZqaesZ9cof7ztW7d2+GDx9Ot27daN26NYMGDeLIkSOnDPd9zTXXFDncd+vWrXG5XIwdO5bY2Ni84b579OhBjRo1TtsXnOG+p06dSps2bZg+fTrPPPMM4DQet27dmlatWtGzZ0/atm3LrFmzaNWqFSkpKWzcuJERI0b49G+Uf7jvtm3b0r17dzZt2sShQ4fo168fbdu25ZJLLjlluO/x48cHpJHaBusrTMYeeHMUbP8aOt8MfR6DiKjAHb+ErdpxgGtfXkrDpDhmjulKQkz4novxT7EG6wvTruA23PepznawPqtBFCb+HBjxHnS9FZa9BNOuhCO/hjqqYmtfryr/+0N7Nv5yhLGvrSQrJ7wb4k0Ju/7DsEsO4Az3nZKSQvv27Rk4cOBZD/edkpJCy5YtOXTokA337deBRFxAvKoe9qFsX+AZIAKYpKqntZSJyO+BBwEF1qjq8DMdM6jDfX//Frw3DmIqw+9fhXpdgvM5JWDOyjT++uYa+rWpxf8NbYfLhuQod2y47/KrRGsQIvKGiFQSkThgPfCDiNxVxD4RwPPA74AWwDARaVGgTGPgXuB8VW0J/NmfOP3WehDc+KnTFfaVfrDsZQjTW3MDOyRz3+XN+PC73Tz0/jobkqOcst97+VOc37m/t5haeGoMVwHzgHrAdUXs0xnYoqpbVTUbmAkMKFDmJuB5VT0AoKp7/IzTfzVbwZjPodHFMO9OePcWOP5bqKMqljE9GzGmZ0Omfbud5z/fEupwTAmLiYkhPT3dkkQ5oqqkp6cTExNzVvv5+xxElIhE4SSI51T1uIgU9VdXB9iZbzkNKHjPpgmAiHyNcxvqQVX9uOCBRGQMMAacpyyDLrYqDJsFXz4BXz4Ov66FIa9B1XOD/9kBdk/fZuw7ksVTCzaRGB/N0M4l8O9nSoXk5GTS0tIKHdLClE0xMTEkJyef1T7+JoiXgG3AGmCRiJwLFNUG4e2md8GkEgk0Bi4CkoGvRKSVqh48ZSfVicBEcNogzjb4YnG54OJ7oXY7eHsMTLwQBk2BRkU/8FKauFzCE4PasP9YNve98z1V4yrQp2XNUIdlSkBUVFTeE8LGnIlft5hU9VlVraOql6tjO3BxEbulAXXzLScDu7yUeU9Vj6vqT8APOAmj9Gja17nlFF8TXhsIXz0ddu0SUREu/ndte9okV+FPM1azdGt6qEMyxpQi/jZS3+5ppBYRmSwiq4CivkovBxqLSAMRqQAMBeYWKPMunkQjIok4t5y2+hNrUFRv5DRet7gKFj4Es6+DzCI7cZUqFStEMnVUJ+pWjeXGV1ew8Zfwit8YEzz+NlKP9jRS9waSgOuBMw7uoqo5wDhgPrABmK2q60TkYRHp7yk2H0gXkfXA58Bdqlo6v95Gxzu3mHo/BhvnwaResHdTqKM6K1XjKvDqDV2IqxDJiMnL2LnfhuQwxvj5HISIfKeqbUTkGeALVX1HRFararvAheiboD4H4aufFsGb10NOFlz9AjQ/ffTI0uyHX44w+MVvSIyP5s2x3ageHx3qkIwxQRbMJ6lXisgC4HJgvogkAO4i9im7GvSEm7+ExMYw6w/w6UNhNXR405oJTBnViZ8P/sboV5ZzNOv0iVmMMeWHvwniBuAeoJOqHgMq4NxmKr8qJzuD/bUfCYufdhqwj+0PdVQ+61i/Gs8Pb8/aXYcZ+9pKsnPKb743przztxeTG6cX0j9E5Cmgu6qG9xjZgRAVA/2fhSufcQb7m3gh7DrzCJWlyaUtavCvq1vz1eZ93PXWGtzu8OqdZYwJDH97MT0O3I4zzMZ64DYR+VcgAisTOoyC6z92bjNN6QOpM0Idkc9+36kud/Vpynupu3j0ww321K0x5ZC/D8pdDqR4ahKIyDRgNc44SgYguQOM+RLeuh7eHQs/r4Q+4yGyQqgjK9ItFzViX0YWU77+iXMqRTP2Qu9zBRtjyqZADPddJd97m9PSm/gkuO5d6DYOlr8M066Aw7tDHVWRRIR/9mtB/7a1efyjjby5YmfROxljygx/E8S/gNUi8oqn9rASGO9/WGVQRKQz6dCgKfDL9067xPZvQx1VkVwu4anBbbmgcSL3vP09CzeE75wYxpiz428j9QygK/C259VNVWcGIrAyq9VAuHEhVIhzahJLJ5b6IToqRLp44Q8daFm7Ere+sYqV28OnV5YxpviKlSBEpH3uC6iFM3bSTqC2Z505kxot4KbP4bxL4aO74J2xkF26n16Oj45kyqhO1Kocy+hXVrD518BMim6MKb2K9SS1iHx+hs2qqiU+tGmpeJL6bLndsOhJ+OJfznwTQ16DqvVDHdUZ7dx/jGte+IZIlzDnj92pXSU21CEZY/xwpiepAzblaKiFZYLItWkBvH0jIDBoslOzKMXW7zrMkJe+pUblGN68uRtV40p/jyxjjHfBnHL0Gi+vXiJyjj/HLXea9HZuOVWqA68NgkVPObWLUqpF7Uq8PLIjO/YfY/S05RzLtiE5jCmLAjHUxiTgWs/rZeAvwNciUtTUoya/6o3gxk+cRuzPHin1Q4d3bVidZ4emsGbnQca9sZrjJ0pvQjPGFI+/CcINNFfVgao6EGgBZOFMIXq3v8GVOxXiYOAk6PMv+OEjePkS2LMx1FEVqm+rWjxyVSs+27iHe+Z8b09bG1PG+Jsg6qtq/o7xe4AmqrofOO7nscsnEeh2C4ycC5kHnfkl1r8X6qgKdW2Xc7nj0ibMWZXG4x+X3mRmjDl7/iaIr0TkAxEZKSIjcWaGWyQiccDBIvY1Z1K/hzNER1IzmD0CPnkATpTOe/239TqP67qey0tfbmXSV6Vv4j9jTPH4OxbTrcA1QA9AgGnAHHXuNRQ1N7UpSuU6cP08+Ohu+Pq/sDsVBk6BuOqhjuwUIsKD/VuSfjSLRz/cQPX4ClzdLjnUYRlj/OTvk9QKLAY+Az4FFqndiA6syGi48r/Q/zlnaI6JF8Ku1d7LTu3nvEIgwiVMGJJCt4bVuevN7/jihz0hicMYEzj+dnP9PbAMGAT8HlgqIoMCEZgpoP11MPpjZ1iOyX1g9euhjug00ZERTBzRgSY1Evjja6tYveNAqEMyxvjB3zaIv+PMJjdSVUcAnYF/+h+W8apOe2dK03pd4L1b4IO/QE52qKM6RUJMFK+M7kRSQjSjX1nOj3szQh2SMaaY/E0QLlXNfy8hPQDHNGcSlwh/eAe63wYrJsMr/eDwrlBHdYpzEmJ4dXRnIlzCiMnL+OVQZqhDMsYUg78X849FZL6IjBKRUcCHwDz/wzJnFBEJvR+Bwa/Ar+vgpQth+zehjuoU9RPjeOX6zhz67Tgjpyzj0DHr9WxMuPG3kfouYCLQBmgLTFRVe0CupLS8Gm5aCDGVYNqVTk2iFPURaFWnMhOv68BP+45y46vLyTx+ItQhGWPOgt+3g1R1jqr+RVXvUNV3AhGUOQvnNIebPoPGveHAVti7ETa8D7+VjsdQup+XyNND2rJi+wHGvbGaHBuSw5iwUdzhvo8A3nYUnN6vlfwN7GyF9WiugeB2w7MpcGgnqBvEBbXbQ8OLoNHFkNzJ6TIbIq9+u43731vHkI51eXxga0QkZLEYY04602iuxXpQTlUT/AvJBJzLBZXrOiPC9vonbP0CfvwcFk+Ar56CqIpw7vlOwmh4EdRo6QzrUUJGdKvPviNZPPvZFpISormzT9MS++xQWze+BwAt71sc4kj8U1bOA+xcfOXvk9TFIiJ9gWeACGCSqj5eYPso4EngZ8+q51R1UokGGa7EBed2d14X3weZh2DbYidhbP0CFvzdKReXdDJZNLzYeWo7yO64rAl7M7J47vMtJMZXYNT5DQotW5b+AxsTrko8QYhIBPA8cBnOVKXLRWSuqq4vUHSWqo4r6fjKnJjK0Kyf8wI4lAZbv4StnzsJ4/s3nfXVGzu3ohpe5IwDFVM54KGICI8MaEV6RjYPfbCeavHR9G9bO+CfY4wJjFDUIDoDW1R1K4CIzAQGAAUThAmGysnQ7lrnpQp71ju3orZ+Aatfg2UTQSKgToeT7Rd1OkJkYGaNi4xw8eywdoyYvIy/zk6lWsUK9GicGJBjG2MCKxQPtdUBduZbTvOsK2igiHwnIm+JSF1vBxKRMSKyQkRW7N27Nxixlm0iTltE93Hwh7fg7u0w6kO44C9OQ/dXT8HU38ET9eH1wfDt/+DX9X53pY2JiuDlkR1plBTPzdNX8H3aocCcjzEmoEJRg/DWMlrwivM+MENVs0RkLM4osZectpPqRJznMOjYsWPpeQAgXEVWcG4v1e8Bl/zD6Sq7bfHJ21GbFzjl4mucbLtoeCFUOvvbRJVjo5g2ujPX/O8bRk1dxlt/7E6DxLhAno0xxk+hSBBpQP4aQTJwylgRqpqeb/Fl4IkSiCv8Xf9hYI8XWwWaX+G8AA7uPNnYvWUhfDfLWZ/Y9GT7xbnnOw/u+aBGpRim39CZQS9+y4gpS5kztjvnVIoJ7DkYY4otFAliOdBYRBrg9FIaCgzPX0BEaqnqbs9if2BDyYZovKpS1xlVtv11znMXe9ad7E67chosfdFpv0julK/9ogNERBV6yIZJ8Uwd1YlhLy9h5NTlzLq5K5ViCi9vjCk5JZ4gVDVHRMYB83G6uU5R1XUi8jCwQlXnAreJSH8gB9gPjCrpOE0RXC6o2dp5df8T5GTBzqUnaxiL/g1fPg4V4p1bVg09NYykpqc9f9G2bhVe/EMHRr+ynJumrWDa6M4hOCFjTEEheQ5CVedRYFA/Vb0/3/t7gXtLOi7jh8hoaNDTefW6H347AD99dbL9YtPHTrmEWvmev7gIEmoC0LNJEv/5fVtun5nKn2emcosKEWLNSsaEUkgShCkHYqtCi/7OC+DA9pO1i80LYM0MZ31S87z2iwHNu7PvihY88sF6JKo3t8bMD1HwxhiwBGFKStVzocNI5+V2w6/fn2y/WDEFlvwPXJHckNyZJue1ZMKPdRh/7CpqT1+Ey+Ui0hVBhEuIiHARGeEiIsJFhMt573JFEOVynbIt0uUiMtLZJypCiHC5iHKJZ9lFhEuIjDj5Psrlylf21DKRLpdTNl8ZG0vKlAfFGqyvNCr3g/WFs+OZsHNJXg1Dd6UiXseCLB63OkdTxPMCxeVl3cn3eFl3+nZxfkrusvNYkYqzHgREcJ9w5sJweRrrC56boJ6O3lr49tPW567D676+HcPb+nzHLbg973P8S46B+N0G8u+jLFilTen+8NfF2jfgg/UZE1BRMSfbJAA5tp+dT19IlB6nZq9xgHoezvP+U1Vxq+J2u/N+qtuN262cUDfqVtzqdtar4nYrbvcJz3tnH3W7Uc19r3nr3Kqo5xjO5zjHcJY97zX/e/cp21Dl+JG9KBAVm3Tq5Vfyvc9b5/xUxLP55HJuudzKS8EyucdQwbMl34VcIO+52LzVrgLX+lNjyP2s3M/57dfNgBJTs9kZf50+pY9A1MD8OMZvu5yBG2Jrt/A/jhD7bdd6clyBGemgIEsQpvSpWI3DLmcsqJrdix6OS3C6w0UEN6piyxt48K63QxyJf/LO448vhDgS/+Wdy83PhzgS/znnEpwZG23+aGOMMV5ZgjDGGOOVJQhjjDFeWYIwxhjjlSUIY4wxXlmCMMYY45UlCGOMMV5ZgjDGGOOVPShnSqWWtSqHOgRjyj1LEKZ0CvTseMaYs2a3mIwxxnhlo7kaY0w5dqbRXK0GYYwxxitLEMYYY7yyBGGMMcarMtMGISJ7ge1+HCIR2BegcEKprJwH2LmURmXlPMDOJde5qprkbUOZSRD+EpEVhTXUhJOych5g51IalZXzADsXX9gtJmOMMV5ZgjDGGOOVJYiTJoY6gAApK+cBdi6lUVk5D7BzKZK1QRhjjPHKahDGGGO8sgRhjDHGK0sQHiLyiIh8JyKpIrJARGqHOqbiEpEnRWSj53zeEZEqoY6puERksIisExG3iIRdl0QR6SsiP4jIFhG5J9TxFJeITBGRPSKyNtSx+EtE6orI5yKywfO3dXuoYyouEYkRkWUissZzLg8F9PjWBuEQkUqqetjz/jaghaqODXFYxSIivYHPVDVHRJ4AUNW7QxxWsYhIc8ANvATcqaphMyKjiEQAm4DLgDRgOTBMVdeHNLBiEJGeQAbwqqq2CnU8/hCRWkAtVV0lIgnASuCqMP29CBCnqhkiEgUsBm5X1SWBOL7VIDxyk4NHHBC2mVNVF6hqjmdxCZAcynj8oaobVPWHUMdRTJ2BLaq6VVWzgZnAgBDHVCyqugjYH+o4AkFVd6vqKs/7I8AGoE5ooyoedWR4FqM8r4BduyxB5CMij4nITuBa4P5QxxMgo4GPQh1EOVUH2JlvOY0wvRCVVSJSH2gHLA1tJMUnIhEikgrsAT5R1YCdS7lKECLyqYis9fIaAKCqf1fVusDrwLjQRntmRZ2Lp8zfgRyc8ym1fDmXMCVe1oVtzbSsEZF4YA7w5wJ3EMKKqp5Q1RScOwWdRSRgtwDL1ZSjqnqpj0XfAD4EHghiOH4p6lxEZCRwBdBLS3lD01n8XsJNGlA333IysCtEsZh8PPfr5wCvq+rboY4nEFT1oIh8AfQFAtKZoFzVIM5ERBrnW+wPbAxVLP4Skb7A3UB/VT0W6njKseVAYxFpICIVgKHA3BDHVO55GnYnAxtU9elQx+MPEUnK7aUoIrHApQTw2mW9mDxEZA7QFKfHzHZgrKr+HNqoikdEtgDRQLpn1ZIw7pF1NfB/QBJwEEhV1T6hjcp3Ih5eLjEAACAASURBVHI58F8gApiiqo+FOKRiEZEZwEU4w0r/CjygqpNDGlQxiUgP4Cvge5z/7wD3qeq80EVVPCLSBpiG8/flAmar6sMBO74lCGOMMd7YLSZjjDFeWYIwxhjjlSUIY4wxXpWZbq6JiYlav379UIdhjDFhZeXKlfsKm5O6zCSI+vXrs2JF2AzTY4wxpYKIbC9sm91iMsYY45UlCGOMMV5ZgihLpvZzXsYYEwBlpg3CGBM4x48fJy0tjczMzFCHYgIkJiaG5ORkoqKifN4nqAnCMybQMziPgU9S1ccLbI8GXgU64AwLMURVt3m2tcGZJKYSzuPwnVTV/lqNKQFpaWkkJCRQv359nKGLTDhTVdLT00lLS6NBgwY+7xe0W0ye2bSeB34HtACGiUiLAsVuAA6o6nnABOAJz76RwGs44yG1xBkD5niwYjXGnCozM5Pq1atbcigjRITq1aufdY0wmG0QvsymNQBnoCmAt4BenpEWewPfqeoaAFVNV9UTQYzVGFPA2SaHIS99y5CXvg1SNMZfxUn2wUwQvsymlVfGM0XmIaA60ARQEZkvIqtE5G/ePkBExojIChFZsXfv3oCfgDHGlGfBTBC+zKZVWJlIoAfO1J89gKtFpNdpBVUnqmpHVe2YlOT1QUBjTJiKj4/Pe9+3b1+qVKnCFVdcccZ9cnJySExM5N577w12eAF36623kpKSQosWLYiNjSUlJYWUlBTeeuutszrOqlWr+PjjjwMSUzAThC+zaeWV8bQ7VMaZGD0N+FJV93kmvJkHtA9irMaYUuyuu+5i+vTpRZZbsGABTZs2Zfbs2QRzKoOcnJyAH/P5558nNTWVefPm0ahRI1JTU0lNTWXQoEFndZxAJohg9mLKm00L+BlnNq3hBcrMBUYC3wKDgM9UVUVkPvA3EakIZAMX4jRiG2NK2EPvr2P9rqKnbF6/2ynjSztEi9qVeODKlj7H0KtXL7744osiy82YMYPbb7+dF154gSVLltCtWzcAli9fzu23387Ro0eJjo5m4cKFVKxYkbvvvpv58+cjItx000386U9/yhu2JzExkRUrVnDnnXfyxRdf8OCDD7Jr1y62bdtGYmIi48eP57rrruPo0aMAPPfcc3Tv3h2Af//730yfPh2Xy8Xvfvc7brrpJgYPHsyqVasA2Lx5M0OHDmXlypU+nf/mzZsZN24c+/btIy4ujkmTJtGkSRNmzpzJo48+SkREBNWqVWPevHk8/PDD/Pbbb3zxxRf84x//OOsEk1/QEoSq5ojIOGA+J2fTWiciDwMrVHUuzrR/0z0zoO3HSSKo6gEReRonySgwT1U/DFasxpjw99tvv7Fw4UJeeuklDh48yIwZM+jWrRvZ2dkMGTKEWbNm0alTJw4fPkxsbCwTJ07kp59+YvXq1URGRrJ///4iP2PlypUsXryY2NhYjh07xieffEJMTAybN29m2LBhrFixgo8++oh3332XpUuXUrFiRfbv30+1atWoXLkyqamppKSkMHXqVEaNGuXzuY0ZM4ZJkybRqFEjvv76a8aNG8eCBQt46KGH+OKLL6hRowYHDx4kNjaW+++/n7Vr1/Lf//7Xj39NR1Cfg/BM4TevwLr7873PBAYXsu9rOF1djTEh5Os3/dyaw6ybuwUznEJ98MEHXHzxxVSsWJGBAwfyyCOPMGHCBH744Qdq1apFp06dAKhUqRIAn376KWPHjiUy0rkMVqtWrcjP6N+/P7GxsYDzMOG4ceNITU0lIiKCTZs25R33+uuvp2LFiqcc98Ybb2Tq1Kk8/fTTzJo1i2XLlvl0XgcPHmTJkiUMHDgwb13uLa7zzz+fESNGMHjwYK655hqfjnc27ElqY0yZMGPGDL7++mtyh/1PT0/n888/55xzzvHaxVNVva6PjIzE7Xamqi743EBcXFze+wkTJlCjRg3WrFmD2+0mJibmjMcdOHAgDz30EJdccgkdOnSgevXqPp2XqpKYmEhqaupp215++WWWLl3KBx98QNu2bfnuu+98OqavbCwmY0zYO3z4MIsXL2bHjh1s27aNbdu28fzzzzNjxgyaNWvGrl27WL58OQBHjhwhJyeH3r178+KLL+Z9G8+9xVS/fv28toE5c+YU+pmHDh2iVq1auFwupk+fzokTzqNavXv3ZsqUKRw7duyU48bExNCnTx/++Mc/cv311/t8blWrVqVWrVq88847ALjdbtasWQPA1q1b6dq1K4888ghVq1bl559/JiEhgSNHjvh8/DOxBGGMKfUuuOACBg8ezMKFC0lOTmb+/PmnbH/77be55JJLiI6Ozls3YMAA5s6di6oya9Ys/vSnP9G2bVsuu+wyMjMzufHGG6lXrx5t2rShbdu2vPHGGwA88MAD3H777VxwwQVEREQUGtMtt9zCtGnT6Nq1K5s2bcqrXfTt25f+/fvTsWNHUlJSeOqpp/L2ufbaaxERevfufVbnP3PmTF588UXatm1Ly5Yt+eCDDwC44447aN26Na1bt+bSSy+lVatWXHLJJaxZs4Z27dqddRfZgiSYXcFKUseOHbXcTxiUO5Lr9daeb/yzYcMGmjdvHuowypynnnqKQ4cO8cgjj4Tk8739XkVkpap29Fbe2iCMMaYEXH311fz444989tlnoQ7FZ5YgjDGmBOS2IYQTa4MwxhjjlSUIY4wxXlmCMMYEhk15W+ZYgjDGGOOVJQhjTKmUO9x3amoq3bp1o2XLlrRp04ZZs2YVuo8N9x0+o7kaY4zfKlasyKuvvkrjxo3ZtWsXHTp0oE+fPlSpUuW0svmH+x4/fnzQpkzNycnJG8MpUJ5//nkAtm3bxhVXXOF1aA1frFq1irVr19K3b1+/Y7IEYYw5s4/ugV++L7rcL55xgHxph6jZGn73uE8f36RJk7z3tWvX5pxzzmHv3r1eE4QN9x0mw32HFXsC2ZiwsGzZMrKzs2nUqNFp22y47xAM9+2Z0+F1VT3g96cZY8KPj9/0g/1Fa/fu3Vx33XVMmzYNl+v05lMb7js0w33XBJaLyCpgCjBfy8oATsaYsHD48GH69evHo48+SteuXb2WseG+QzDct6r+A2iMM/vbKGCziIwXkdPreMYYE2DZ2dlcffXVed+UvbHhvkM43LenxvCL55UDVAXeEpF/ByQKExDrdh9i3e5DoQ7DmICaPXs2ixYt4pVXXsnr+lnw27QN9x2i4b5F5DZgJLAPmAS8q6rHRcQFbFbVUlGT8Gu47zLSSL1ufA8AWt63OMSRmHBnw30HR1kc7jsRuEZVt+dfqapuEbniTDuKSF/gGSACmKSqjxfYHg28CnQA0oEhqrot3/Z6wHrgQVV9CmOMCVNldbjveUBe/y8RSQBaqOpSVd1Q2E4iEgE8D1wGpOE0dM9V1fX5it0AHFDV80RkKPAEMCTf9gnARz6fjTHGlFJldbjvF4CMfMtHPeuK0hnYoqpbVTUbmAkMKFBmADDN8/4toJd4mv9F5CpgK7DOh88yxgSYdVYsW4rz+/QlQUj+bq2q6sa3mkcdYGe+5TTPOq9lVDUHOARUF5E44G7goTMGJjJGRFaIyIq9e/f6EJIxxhcxMTGkp6dbkigjVJX09PS8rri+8uVCv9XTUJ1ba7gF55t9UbwNglLwr62wMg8BE1Q140xjqajqRGAiOI3UPsRkTMkLw04QycnJpKWlkf+L1/FDvwAQVblmqMIKmPJ4LjExMSQnJ5/VsX1JEGOBZ4F/4Fy8FwJjfNgvDaibbzkZ2FVImTQRiQQq47R3dAEGebrRVgHcIpKpqs/58LnGGD9FRUXRoEGDU9atG38TAM3LQC85OxffFJkgVHUPMLQYx14ONBaRBsDPnmMML1BmLk4X2m+BQcBnnttZF+QWEJEHgQxLDsYYU7J8GYspBqe3UUsg7waWqo4+036qmuMZx2k+TjfXKaq6TkQeBlao6lycp7Oni8gWnJpDcRKRMcaYIPDlFtN0YCPQB3gYuBYotHtrfqo6D6ebbP519+d7nwl4f3b+ZJkHffksU8aE4X17Y8oaX3oxnaeq/wSOquo0oB/QOrhhGVN22BAoJlz5kiCOe34eFJFWOA3J9YMWkTHGmFLBl1tME0WkKk4vprlAPPDPoEZljDEm5M6YIDwD8h32TBa0CGhYIlEZY4wJuTPeYvI8NT2uhGIxxhhTivjSBvGJiNwpInVFpFruK+iRmbNz4jix7mPEuw/D8cyiyxtjTBF8aYPIfd7h1nzrFLvdFFqZhyFtOexYAju+hZ9X0jDHmcGK/zSBltdAynBI7gRnGK7EGGMK48uT1A2KKmNKwKGfYeeSkwnh13WgbhAX1GwD7Uewc+U83Lg4t0kX+G4WrJwK1RpB22HQdghUqRfqszDGhBFfnqQe4W29qr4a+HAMAG437N3oJIIdnqRwaIezLSoOkjtCz79Bva7O++gEAA6v/topc81EyDoC6+fCmhnw+aPOq/4FTrJo0T9vH2OMKYwvt5g65XsfA/QCVuHMBGcC4Xgm7Fp1MiHsXAqZnger4ms4iaDbLc7PGq0hwodfW3QCtLvWeR3Y7tQo1syA926BeXdC8yudZNGgJ7gKn3c3VHIfLGsZ4jiMKc98ucX0p/zLIlIZZ/gNU1zH9nsSgad2sGs1nMh2tiU2hRZXQb1uUK8LVG3gfxtC1XPhwr9Bz7tg5zInUax920kalepAm99D2+GQ1MT/czPGlBm+1CAKOgY0DnQgZZYqHPgJdiw9WUPY94OzzRUFddpD1z9C3a5QtwvEVQ9eLCJO0qnXBfo+Dj/MgzUz4etnYfEEqN3eadhuNRAqWkc1Y8o7X9og3ufkRD8uoAUwO5hBhbUTOfDr9ycbk3csgYxfnW0xlZ0k0HaIU0Oo3Q6iYkMTZ1QMtLrGeR35Fb5/06lZzLsTPr4XmvRxksV5l0FkhdDEaIwJKV9qEE/le58DbFfVtCDFE36yMk7tbpq2Ao4fdbZVrgcNLnTaDup1g6Rm4PLl0ZMSllADuo9zXr98D6kz4PvZsPEDqFgdWg+GtkOhVop1mTWmHPElQewAdnuG5kZEYkWkvqpuC2pkpdXh3ad2N/1lLegJQKBmK6dRuF5X55ZR5YJTcIeBmq2hb2u47CH48TNIfQNWTIGlL0JScydRtBkClWqFOlJjTJD5kiDeBLrnWz7hWdfJe/EyxO2GfZucRLDT04ZwYJuzLTLW6WJ6wV+de/rJnZxbSGVFRJRzm6lJH/jtAKx7x6lZfPoALHwIGl7s9IJq1g8qVAx1tMaYIPAlQUSqanbugqpmi0jZvCmdk+X0KMp99mDnEufiCBCX5NQMOt3k3C6q1ca5iJYHsVWh42jntW8LfDfTadx++0aokAAtr3KSxbnd7RaUMWWILwlir4j090wRiogMAPYFN6wSdHCHUyvIOgz/qgsnspz11c9zvh3X6+a8qjW0ix9A4nlwyT/govtg+2InUax9G1ZPhyrnep7aHgrV7AF8Y8KdLwliLPC6iDznWU4DvD5dHZayj8Lhn6FCPHT21A7qdoH4pFBHVrq5XM5Ddg16wuVPwob3nV5QXz4BXz7u/Du2HebULsrSrTdjyhFfHpT7EegqIvGAqOqR4IdVghKbOg3Krgjo81ioo/HLPncCG0/UJurXIzRMjCMyooR6TFWIc2oNbYfCoTTnAbzUGfD+bfDR35yaWNvh0PAi354CN8aUCr48BzEe+LeqHvQsVwX+qqr/8GHfvsAzQAQwSVUfL7A9GmfIjg5AOjBEVbeJyGXA40AFIBu4S1U/O6sz85XLVSqHmjgbRzKP8+KXP/Jyxs1kE8W/JiwiJspFs5qVaFWnEi1rV6ZV7co0qRlPdGSQz7VystNw3+Mv8PMqWPMGfP8WrJ0D8TWhzWAnWdRoEdw4jDF+8+Xr3O9U9b7cBVU9ICKX40xBWigRiQCeBy7DuS21XETmqur6fMVuAA6o6nkiMhR4AhiC08Zxparu8syDPR8Iwz6jwXX8hJsZy3bwzKebST+azUWRm+hXYSWufv9h3a7DrN11iPdW7+K1Jc5Af5EuoXGNBFrVrkTL2pVoVacyzWtVIi46CN/qRSC5g/PqMx42zXduQS15Ab75P6jV1rkF1XowxCUG/vONMX7z5coQISLRqpoFznMQQLQP+3UGtqjqVs9+M4EBQP4EMQB40PP+LeA5ERFVXZ2vzDogJn8M5Z2qMn/dr/z7441s3XeUrg2rMfXy5kS8+i8AWnZIZmAHp6zbrezYfywvYazbdZjPNu7hzZXOs44i0CAxjla1K+cljZa1K1GlYgA7qkVGOyPItugPGXud2sSaN+Dje2DBP6Bxb+f2VJO+TtnSQNXp1Xb8GGRnQPYx5wHI7GMF1h1z2rGOH3OWszNOvveUb3R8M6Dwv26hPisA3KqowglV3G7Frc46d4HlE27POs9yxexMVIWfHm4T6lPwW2yO0zGzrJzLdqkVlIEtfUkQrwELRWSqZ/l6YJoP+9UBduZbTgO6FFZGVXNE5BBQnVN7SQ0EVntLDiIyBhgDUK9e+ZjrYPWOA4yft4Hl2w7QKCmOSSM60qv5OYgI67yUd7mE+olx1E+Mo18b5+E2VeXXw1ms/flQXtJYsW0/c9fsytuvTpXYvISRe5vqnIRoxN+eXPFJ0HWs8/p1vVOr+G62My5UTBVnHKiU4c4FuqjPUnUGOSziAu1sP3rqBT1vn6MFyh89+V7dZ3Fi4rTFRFV0fua9r0g2TrKNqd6o2P9suRf0nBNKjlvJcbs9791560643RzPfZ+7za2nlXOrFv2BODXOyAiX89MloIdw4SYqNvxrfDlH9gIQGRv+nVFyjuwlp1jD6hXNl0bqf4vId8ClgAAfA+f6cGxv/7sL/mWesYyItMS57dS7kNgmAhMBOnbs6NtffZjakX6MJ+Zv5MPvdpMYH81jV7diSMe6xWqIFhFqVo6hZuUYLm1RI2/9/qPZrPMkjLU/H2L9rsMsWP9r3vbE+GhP0jjZrlG3Wmzxk0aNFtD7Eej1AGz9wkkWqa/DismcRzS/uWJh1h9OvWgXvMjribP7zKg458G+3At57s+K1fMu6E6ZuHzvc8vHF74uMsZrQnO7lfX/uowMjaHWBc+TkZnD0ewcjmTmkJGVQ4bn55HMHI5medZlndx+NLdMdg6+XNcrRLpIiI4kPsbziva8Yk7+TIiOJM6zPiEmkvjoqNPKVoyKwOU69XzWje8BuGh51wdn929eCjnnQpk5lzh+C8qxfU07vwBu4PfAT8AcH/ZJA+rmW04GdhVSJk1EIoHKwH4AEUkG3gFGeHpSlUsHjmbzf59tYfqSbUS6XNx2yXmMubAR8UFoN6gWV4ELGidxQeOT36oysnLYsNtJGGt/Psy6XYdYvGUfJ9zO1SohJtJJGrUr07KO87NhUjwRrrNIGhGR0PhS55V5CNa9S84H91LRfdR5MC/3Yhxbx+s39JMXdG/r8l3QI2MDMhZWzgk3+49ls+9INvv2ZrEvI530jGz2ZWSxNyOLfRnZ7DuSRfrRLNIzsslx3+Hs+Oxir8cTwblYey7OcZ4Ld+0qMZ6LdhTx0RGei3hUIRd5Z7lCZCkc68uErUKvMiLSBBgKDMPpYTQLp5vrxT4eeznQWEQaAD97jjW8QJm5wEjgW2AQ8JmqqohUAT4E7lXVr8/ifMqMzOMnePXbbTz32RYysnIY3KEuf+ndhBqVYko0jvjoSDrVr0an+ieH/848foJNvx7JSxhrdx1m+pLtZOU4t2Riolw0r1UpL3G0qlOZxjV87EEVUxk6jGTb/JcBaHmr94tqoGXlnMi7yO/LyHIu/kc9P3PXeS7+B45le/02XyHSRVJ8NInxFahZOYZWdSqRGB/N8W9fIl4yaTrogbyLef5v7BUrRPh/686YIDjT19CNwFc4vYm2AIjIHb4e2NOmMA6nB1IEMEVV14nIw8AKz5PZk4HpIrIFp+Yw1LP7OOA84J8i8k/Put6quucszi0sud3K+9/t4t8f/8DPB3/jwiZJ3Ht5M5rVrBTq0PLEREXQJrkKbZKr5K3LOeHmx71HnYThSRwFe1A1qZFwSkN40HpQeRzNyiE9I9vzrf7khT/9aL4k4PnWfyQzx+sx4ipEkJgQTWJ8NPWrx9GxfjUS46NJiq9AYnw0iQnRVI+rQGJCNAnRkV4v9OtWLQegZauaQTtXY4LhTP87B+JcsD8XkY+BmXhvMyiUqs4D5hVYd3++95nAYC/7PQo8ejafVRYs2ZrO+Hkb+C7tEC1qVeKJgW3o0Tg8GgQjI1w0rZlA05oJXNPeWXc2Pahy2zXO1INKVTn8W47nm73nVk6+b/b7CiSC3457b5+oHBtFoucC37xWJS7Id7FP9NQAnJ/RxFYI72dkjPFHoQlCVd8B3hGROOAq4A6ghoi8ALyjqgtKKMYyb8ueIzz+0UY+3bCHWpVj+M/gtlzdrs5pjYThpqgeVLmJw1sPqnrHriFOsjgxdVnexT89I5vsE6f3LHKJ036Se1E/t15Fqsfnu9gnRJMUH031+ApUj4u2+/TG+MiXXkxHgddxxmOqhvON/x7AEoSf9h7J4r+fbmLm8p3ERkVwV5+m3NCjATFRZfdbq689qFavTSRTo6iVkUVifDRNayZQPb6C5x6/55XgXPCrxVU4u0ZxY4xPzuoGsKruB17yvEwxHcvOYdJXP/HSlz+SlePmD13qcVuvxlSPLyUPiYVAwR5U68bfBkDLP5VMI7Ux5nQ2cloJOuFW3lq5k6c/2cSvh7Po07IGd/dtRsOk+FCHZowxp7EEUQJUlS837eXxjzay8ZcjpNStwnPD25/SddQYY0obSxBBtm7XIf41byOLt+yjXrWKPD+8PZe3rmn93o0xpZ4liCDZfeg3npq/ibdXp1E5Nop/XtGCP3StF/zhto0xJkAsQQTYkczjvPDFj0xe/BMKjLmgIbdcfB6VY8vJ/NXGmDLDEkSA5M7N8N9PN7P/aDZXpdTmzj5NSa5aMdShGWNMsViC8JO3uRnuu7z5KcNQGGNMOLIE4Yf8czOcd048k0d25JJm51gDtDGmTLAEUQyBnJvBGGNKK0sQZ+G0uRl6NWZMz4ZBmZvBGGNCza5sPig4N8PvO9bljstKfm4GY4wpSZYgzqDg3AwXNU3i3t81p2nNhFCHZowxQWcJohAF52b496A2nH9eeMzNYIwxgWAJAli3+xAALSm7czMYY8zZsgThccBdkb+/8z0zl++kYlQEf+vblNHnl+25GYwx5kzKfYI4lp3DjKzzmZPVhePLd9rcDMYY4xHUjvsi0ldEfhCRLSJyj5ft0SIyy7N9qYjUz7ftXs/6H0SkT7BiXL/rMK9l9aRd5E8suKMnDw1oZcnBGGMIYg1CRCKA54HLgDRguYjMVdX1+YrdABxQ1fNEZCjwBDBERFoAQ3GaBWoDn4pIE1X1Pgu9HzrWr8YLcROpF5FOw6S/BvrwxhgTtoJZg+gMbFHVraqaDcwEBhQoMwCY5nn/FtBLnHEqBgAzVTVLVX8CtniOFxT1ItKDdWhjjAlboqrBObDIIKCvqt7oWb4O6KKq4/KVWespk+ZZ/hHoAjwILFHV1zzrJwMfqepbBT5jDDAGoF69eh22b98elHMxxpiySkRWqmpHb9uCWYPw1i+0YDYqrIwv+6KqE1W1o6p2TEpKKkaIxhhjChPMBJEG1M23nAzsKqyMiEQClYH9Pu5rjDEmiIKZIJYDjUWkgYhUwGl0nlugzFxgpOf9IOAzde55zQWGeno5NQAaA8uCGKsxxpgCgtaLSVVzRGQcMB+IAKao6joReRhYoapzgcnAdBHZglNzGOrZd52IzAbWAznArcHowWSMMaZwQWukLmkishfwp5U6EdgXoHBCqaycB9i5lEZl5TzAziXXuarqtRG3zCQIf4nIisJa8sNJWTkPsHMpjcrKeYCdiy9sCjRjjDFeWYIwxhjjlSWIkyaGOoAAKSvnAXYupVFZOQ+wcymStUEYY4zxymoQxhhjvLIEYYwxxitLEB4i8oiIfCciqSKyQERqhzqm4hKRJ0Vko+d83hGRKqGOqbhEZLCIrBMRt4iEXZfEouZECRciMkVE9ngG2AxrIlJXRD4XkQ2ev63bQx1TcYlIjIgsE5E1nnN5KKDHtzYIh4hUUtXDnve3AS1UdWyIwyoWEemNM2xJjog8AaCqd4c4rGIRkeaAG3gJuFNVV4Q4JJ955kTZRL45UYBhBeZECQsi0hPIAF5V1VahjscfIlILqKWqq0QkAVgJXBWmvxcB4lQ1Q0SigMXA7aq6JBDHtxqER25y8IjDy+ix4UJVF6hqjmdxCc5gh2FJVTeo6g+hjqOYfJkTJSyo6iKc4XDCnqruVtVVnvdHgA1AndBGVTzqyPAsRnleAbt2WYLIR0QeE5GdwLXA/aGOJ0BGAx+FOohyqg6wM99yGmF6ISqrPNMctwOWhjaS4hORCBFJBfYAn6hqwM6lXCUIEflURNZ6eQ0AUNW/q2pd4HVg3JmPFlpFnYunzN9xBjt8PXSRFs2XcwlTPs1rYkJDROKBOcCfC9xBCCuqekJVU3DuFHQWkYDdAgzaaK6lkape6mPRN4APgQeCGI5fijoXERkJXAH00lLe0HQWv5dwY/OalFKe+/VzgNdV9e1QxxMIqnpQRL4A+gIB6UxQrmoQZyIijfMt9gc2hioWf4lIX+BuoL+qHgt1POWYL3OimBLmadidDGxQ1adDHY8/RCQpt5eiiMQClxLAa5f1YvIQkTlAU5weM9uBsar6c2ijKh7P/BrRQLpn1ZIw7pF1NfB/QBJwEEhV1T6hjcp3InI58F9OzonyWIhDKhYRmQFchDOs9K/AA6o6OaRBFZOI9AC+Ar7H+f8OcJ+qzgtdVMUjIm2AaTh/Xy5gtqo+HLDjW4Iwxhjjjd1iMsYY45UlCGOMMV5ZgjDGGOOVJQhjjDFeWYIwxhjjlSUIY86CiGQUXeqM+78l9mlsswAAAhNJREFUIg097+NF5CUR+dEzEuciEekiIhU878vVg6ym9LEEYUwJEZGWQISqbvWsmoQzAF5jVW0JjAISPQP7LQSGhCRQYzwsQRhTDOJ40jNm1PciMsSz3iUi//PUCD4QkXkiMsiz27XAe55yjYAuwD9U1Q3gGfX1Q0/Zdz3ljQkZq8IaUzzXAClAW5yni5eLyCLgfKA+0Bo4B2co6Smefc4HZnjet8R5KvxEIcdfC3QKSuTG+MhqEMYUTw9ghmckzV+BL3Eu6D2AN1XVraq/AJ/n26cWsNeXg3sSR7ZnQpv/b+/uURoIwjCO/99GLeyC6AG09gCC4BksvIBnEKxEPIZn8KPQNoKdYuMXNtpYCVZWESzGYmYRZZJiWQwJ/1+VDeyyW4SHySzvI42FASG1UxvlPep7gAEwVz4/AqsRMeo3OAt8trg3qRMGhNTOJbBVyloWgHXgmlz5uFn2IhbJA+4aT8AyQErpBbgB9st0USJipenAiIge8J5S+vqvB5L+MiCkdk6AO+AW6AM75S+lI3IPxAO5R/sK+CjnnPM7MLaBJeA5Iu6BQ376IjaAiZsuquniNFepYxExX0rke+RVxVpK6a3M678ox8M2p5trHAO7E9zHrSngW0xS985KicsMcFBWFqSUBhGxR+6lfh12cikXOjUcNG6uICRJVe5BSJKqDAhJUpUBIUmqMiAkSVUGhCSp6hvO1ioyR/MQrwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#pd.DataFrame(grid.cv_results_).to_csv('LogisticGridSearchCV_Otto.csv')\n",
    "#cvresult = pd.DataFrame.from_csv('LogisticGridSearchCV_Otto.csv')\n",
    "#test_means = cv_results['mean_test_score']\n",
    "#test_stds = cv_results['std_test_score'] \n",
    "#train_means = cvresult['mean_train_score']\n",
    "#train_stds = cvresult['std_train_score'] \n",
    "\n",
    "\n",
    "# plot CV误差曲线\n",
    "test_means_neglogloss = grid_neg_log_loss.cv_results_[ 'mean_test_score' ]\n",
    "test_stds_neglogloss = grid_neg_log_loss.cv_results_[ 'std_test_score' ]\n",
    "\n",
    "test_means_accuracy = grid_accuracy.cv_results_['mean_test_score']\n",
    "test_means_accuracy = grid_accuracy.cv_results_['std_test_score']\n",
    "\n",
    "# plot results\n",
    "n_Cs = len(Cs)\n",
    "number_penaltys = len(penaltys)\n",
    "test_scores_neglogloss = np.array(test_means_neglogloss).reshape(n_Cs,number_penaltys)\n",
    "test_stds_neglogloss = np.array(test_stds_neglogloss).reshape(n_Cs,number_penaltys)\n",
    "test_scores_accuracy = np.array(test_means_accuracy).reshape(n_Cs, number_penaltys)\n",
    "test_stds_accuracy = np.array(test_means_accuracy).reshape(n_Cs, number_penaltys)\n",
    "\n",
    "x_axis = np.log10(Cs)\n",
    "plt.subplot(2,1,1)\n",
    "for i, value in enumerate(penaltys):\n",
    "    #pyplot.plot(log(Cs), test_scores[i], label= 'penalty:'   + str(value))\n",
    "    plt.errorbar(x_axis, -test_scores_neglogloss[:,i], yerr=test_stds_neglogloss[:,i] ,label = penaltys[i] +' Neglogloss Test')\n",
    "    #plt.errorbar(x_axis, -train_scores[:,i], yerr=train_stds[:,i] ,label = penaltys[i] +' Train')\n",
    "plt.legend()\n",
    "plt.xlabel( 'log(C)' )                                                                                                      \n",
    "plt.ylabel( 'logloss' )\n",
    "plt.subplot(2,1,2)\n",
    "for i, value in enumerate(penaltys):\n",
    "    plt.errorbar(x_axis, test_scores_accuracy[:, i], yerr=test_stds_accuracy[:,i], label = penaltys[i] + ' Accuracy Test')\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel( 'log(C)' )                                                                                                      \n",
    "plt.ylabel( 'Accuracy' )\n",
    "plt.savefig('LogisticGridSearchCV_C.png' )\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
